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I n t roduc t ion 

In design of engineering systems, the "what if" questions often arise 
such as: what will be the change of the aircraft payload, if the wing 
aspect ratio is incremented by 10 % ? Answers to such questions are 
commonly sought by incrementing the pertinent variable, and reevaluating 
the major disciplinary analyses involved. These analyses are contributed 
by engineering disciplines that are, usually, coupled, as are the 
aerodynamics, structures, and performance in the context of the question 
above . 

The "what if" questions may be answered exactly without using the 
" increment-and-reevaluate" approach and without finite differencing of 
the system analysis, by using a method introduced in ref.l for 
calculating the first derivatives of behavior of coupled systems with 
respect to design variables. The method called the system sensitivity 
analysis has been shown applied in formal optimization that used the 
derivatives to guide search in multidimensional design space, e.g., 
ref. 2 and 3. Obviously, if the problem is strongly nonlinear, the 
efficiency of such search will improve if second and, possibly, higher 
order derivatives are available to the search algorithm. A need for the 
second derivatives to enhance the method of ref.l applied in such 
problems is discussed in ref. 4. 

This note extends the algorithm of ref.l to include the derivatives of 
the second and higher orders, again, without finite differencing of the 
system analysis. It achieves that by recursive application of the same 
implicit function theorem that underlies ref.l. As a supplement to 
ref.l, the note is not self-contained, it references equations in ref.l 
(such references are shown as "eq.no/ref . 1" ) and requires that reference 
as a prerequisite. 

First Order Sensitivity Analysis 

Sensitivity problem stated in ref.l calls for calculation of the 
derivatives of a vector Y solving the governing equations, eq. 1/ref. 1, 
with respect to a design variable Xk. The algorithm developed in ref.l 
yields the derivative of Y as a solution vector Z of the sensitivity 
equations which are linear, simultaneous, algebraic equations of the 
form 

1 A Z = R 

In ref.l, the terms in the above equation are defined by two equivalent 
sets of equations: either eq. 4/ref. 1 (based on the residuals and called 
the Global Sensitivity Equations 1, GSE1) or eq. 8/ref. 1 (based on the 
output/input partial sensitivity derivatives and called the Global 
Sensitivity Equations 2, GSE2). The contents of the matrix of 
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coefficients A and the right hand side (RHS) vector R are different for 
the above two alternative formulations but the solution vector Z is the 
same and represents the derivative of Y with respect to the k-th element 
of the vector of design variables X. It has the meaning of the total 
derivative because it reflects both the direct and indirect influences 
of Xk on Y. 

Sensitivity Analysis of Second and Higher Orders 

Generalization of the above first order sensitivity analysis to the 
higher orders is straightforward by taking advantage of the linearity of 
eq.l above. Although most of the linear algebra texts ignore the 
matter, algorithms for sensitivity analysis of the linear algebraic 
equations solution have been developed in structural sensitivity 
analysis. They stem from the same implicit function theorem that was the 
basis for ref.l, e.g., ref. 5 and 6, and they extend to the second order 
derivatives, e.g., ref. 7 and 8. The pattern established in these 
algorithms will be adapted to solve the problem at hand. 

A compact notation for the derivatives will be used in the remainder of 
the note 

3 ()<Jklm...= 3<i( )/ 9Xk 3X 1 3Xm . . . 

where any subscript may by repeated as required to form a high order, 
mixed derivative with respect to any combination of variables X. 

In the above notation, the correspondence of eq.l above to eq . 4/ref. 1 
and eq. 8/ref. 1 makes the derivatives of Z with respect to X equivalent 
to the derivatives of Y with respect to X, as follows (Z is already the 
first derivative of Y) : 


ZO 

= Yik 

Zii 

= Y 2 kl 

Z 2 1 m 

= Y3klm 


ZNlm. . . = YN+lklm. . . 

Repeated differentiation of eq.l yields the derivatives of Z according 
to a regular pattern, shown below up to the fourth derivative 
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Z 2 1 m 

= R 2 1 m 

- Aim 

Zll 

- D 1 m ( A 1 1 

ZO) ; 


7 

A 

Z 3 lmn 

= R 3 lmn 

1 

> 
►— * 
3 

Z 2 1 m 

- D 1 n ( A 1 m 

Zll) - 

D 2 mn (All ZO); 
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A 

Z* Imnp 

= R*lmnp 

- Aip 

Z 3 lmn 

- D 1 p ( A 1 n 

Z 2 1 m ) 

- D 2 np(Alm Zll) 


- D3mnp(All ZO) ; 


etc . 
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where DsimnC ) is a shorthand for the q-th mixed derivative of the 
product of the pair of functions named in the parentheses, obtained by 
the usual rules of the product differentiation. Here, again, any 
subscript may be repeated. Once the derivatives of Z are obtained, the 
derivatives of Y are available from eq.4. 

The computational cost of a full set of mixed, derivatives of the N-th 
order obtained by means of the above algorithm escalates rapidly with 
the value of N. Firstly, because the number of derivatives in a complete 
set is proportional to the number of the design variables d raised to 
the N-th power, and, secondly, because the eq.5-8 are recursive. 
Inspection of the A terms and the D terms reveals that in the 
calculation of a mixed derivative of the N-th order referencing a 
particular variable, the recursivity requires derivatives of the orders 
N-l and lower with respect to the same variable as prerequisites. For 
example, to solve eq.8 for Z 4 lmnp, the prerequisites are the derivatives 
Z 3 lmn, Z 3 lmp, Z 3 lnp, and Z 3 mnp. 

On the other hand, due to the recursivity all the eq.5-8 share the same 
matrix A so that matrix needs to be factored only once and reused when 
the equations are solved one after another. Each solution requires, 
then, only a backsubstitution of its RHS over the factored A - a task 
ideally suited to the modern vector and parallel processor computers. 

Special Case of Single Variable 

The recursivity becomes also an advantage if there is only one variable 
with respect to which the derivatives are to be taken. It is the case of 
interest because it occurs in searching design space along a line 
defined by a search direction S - an operation that is a part of many 
optimization algorithms. Accurate extrapolation based on higher order 
derivatives may reduce the need for costly repetitions of the system 
analysis in that operation by widening the move limits. Assuming that 
the direction vector S has already been generated, all the variables Xk 
become linked to the step length h through S so that 

9 Xnew — X old + h S 

This relation enables one to substitute X with h in eq . 1/ref. 1, and by 
chain differentiation to replace the right hand sides of eq. 4/ref. 1 and 
eq . 8/ref . 1 with 

d 

10 RHS(h) = Z. Sk RHS(Xk) 

fce*l 

where RHS(h) is the right hand side vector reflecting the presence of 
only a single variable h, RHS(Xk) refers to the RHS in eq. 4/ref. 1 or 
8/ref. 1 related to a variable Xk, and Sk is the k-th element of S. 

With the problem reduced to a single variable h, the variable 
identification subscripts are no longer needed, and R in eq.l may be 
replaced with 

11 R = RHS(h) 

from eq.10 to transform the pattern of equations, eq . 5 through eq.8, to 
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A 

Zi = 

Ri - A 1 

Z°; 
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A 

Z2 = 

R2 - A 1 

Zi - 

Di(Ai 

Z°) ; 
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A 

Z3 = 

R3 - A 1 

Z 2 - 

D 1 ( A 1 

Z 1 ) - D 2 ( A 1 Z0); 
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A 

Z< = 

R< - A 1 

Z 3 - 

D 1 ( A 1 

Z 2 ) - D 2 ( A 1 Zi) 
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A 

ZN = 

RN - Ai 

ZN- 1 

-ii 

s-y 

D 3 ( A 1 ZN-i-s-i) 

By the 

Leibniz ' s 

rule (e. 

g. . 

c 

ref .9) 

the D-terms may 
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D s ( A 

.1 ZN- 

■ X -s - 1 ) = 

£ 

4 s O 

( C 3 q A 1 + a_< J ZN-s-l+q) 

where 

Cib are the binomial coefficients 
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C 3 q 

= e ! 

/ (q! (s 

-q) ! 

) ; C 3 a 

= C 3 o = 1; C 3 i = 


also obtainable from the Pascal's triangle (ref. 9)'. 


Structural Sensitivity Analysis: A Special Case of A Linear 

System 

Another special case arises when the governing equations, eq . 1/ref. 1, 
are linear algebraic simultaneous equations. For an instance, consider 
the load-deflection equations in the displacement method of structural 
analysis 

19 K W = P 

where the stiffness matrix K and the load vector P may be functions of 
design variables X, so that the vector of displacements becomes an 
implicit function of X. If one substitutes the terms in eq . 1 as follows 

20 A = K ; Z - 3W/ 3Xk; R = 9P/ 3Xk - 3K/ 3Xk W; 

then eq.4 through eq.18 formulate the structural sensitivity analysis 
through the N-th order derivatives (including the case of derivatives 
taken in a direction S in the design space). 

Since K is a linear function for a broad class of design variables (e.g, 
the wall thickness in a membrane structure, or the moment of inertia of 
a beam cross-section), the patterns of eq.5 - 8 and 12 - 16 simplify 
drastically because the terms comprising the second and higher 
derivatives of K vanish. When one wishes to compute derivatives for 
only a subset of W, then further simplifications and computational cost 
savings are available by the use of an alternative known as the adjoint 
variable method (ref. 6 and 7 discuss that method including the second 
derivative applications) that is, however, beyond the scope of this 
note . 
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